#pragma once

#include <Eigen/Core>

namespace scipy {
namespace ndimage {
typedef enum {
    NI_EXTEND_FIRST         = 0,
    NI_EXTEND_NEAREST       = 0,
    NI_EXTEND_WRAP          = 1,
    NI_EXTEND_REFLECT       = 2,
    NI_EXTEND_MIRROR        = 3,
    NI_EXTEND_CONSTANT      = 4,
    NI_EXTEND_GRID_WRAP     = 5,
    NI_EXTEND_GRID_CONSTANT = 6,
    NI_EXTEND_LAST          = NI_EXTEND_GRID_WRAP,
    NI_EXTEND_DEFAULT       = NI_EXTEND_MIRROR
} NI_ExtendMode;

NI_ExtendMode extend_mode_to_code(std::string mode)
{
    // Convert an extension mode to the corresponding integer code.

    if (mode == "nearest")
    {
        return NI_EXTEND_NEAREST;
    }
    else if (mode == "wrap")
    {
        return NI_EXTEND_WRAP;
    }
    else if (mode == "reflect")
    {
        return NI_EXTEND_REFLECT;
    }
    else if (mode == "grid-mirror")
    {
        return NI_EXTEND_REFLECT;
    }
    else if (mode == "mirror")
    {
        return NI_EXTEND_MIRROR;
    }
    else if (mode == "constant")
    {
        return NI_EXTEND_CONSTANT;
    }
    else if (mode == "grid-wrap")
    {
        return NI_EXTEND_GRID_WRAP;
    }
    else if (mode == "grid-constant")
    {
        return NI_EXTEND_GRID_CONSTANT;
    }
    else
    { // raise RuntimeError('boundary mode not supported')
        return NI_EXTEND_REFLECT;
    }
}
/// @brief Calculate a 1-D correlation along the given axis.
/// The lines of the array along the given axis are correlated with the
/// given weights.
/// @param input
/// @param weights
/// @param mode
/// @param cval
/// @param origin
/// @return
Eigen::RowVectorXd correlate1d(Eigen::RowVectorXd input,
                               Eigen::RowVectorXd weights,
                               std::string mode = "reflect",
                               float cval       = 0.0,
                               int origin       = 0)
{
    eigen_assert((weights.size()) && "no filter weights given");
    if (!weights.size())
    {
        return {};
    }
    int weights_size = weights.size();
    eigen_assert(((-weights_size / 2) <= origin) && (origin <= ((weights_size - 1) / 2)) &&
                 "Invalid origin; origin must satisfy "
                 "-(len(weights) // 2) <= origin <= "
                 "(len(weights)-1) // 2");
    if (!(((-weights_size / 2) <= origin) && (origin <= ((weights_size - 1) / 2))))
    {
        return {};
    }
    int symmetric = 0;
    /* test for symmetry or anti-symmetry: */
    const int filter_size = weights.size();
    const int size1       = filter_size / 2;
    const int size2       = filter_size - size1 - 1;
    double *fw            = weights.data();
    if (filter_size & 0x1)
    {
        symmetric = 1;
        for (int ii = 1; ii <= (filter_size / 2); ++ii)
        {
            if (std::fabs(fw[ii + size1] - fw[size1 - ii]) > __DBL_EPSILON__)
            {
                symmetric = 0;
                break;
            }
        }
        if (symmetric == 0)
        {
            symmetric = -1;
            for (int ii = 1; ii <= (filter_size / 2); ++ii)
            {
                if (std::fabs(fw[size1 + ii] + fw[size1 - ii]) > __DBL_EPSILON__)
                {
                    symmetric = 0;
                    break;
                }
            }
        }
    }
    const int size1_extern = size1 + origin;
    const int size2_extern = size2 - origin;
    const int line_length  = input.size();
    Eigen::RowVectorXd input_extend =
        Eigen::RowVectorXd::Zero(line_length + size1_extern + size2_extern);
    for (int cc = 0; cc < input.size(); ++cc)
    {
        input_extend[cc + size1_extern] = input[cc];
    }
    if (size1_extern + size2_extern > 0)
    {
        NI_ExtendMode extend_mode = extend_mode_to_code(mode);

        int size_before = size1_extern;
        int size_after  = size2_extern;
        double *buffer  = input_extend.data();
        double *first   = buffer + size1_extern;
        double *last    = first + line_length;
        double *src, *dst, val;

        if ((line_length == 1) && (extend_mode == NI_EXTEND_MIRROR))
        {
            extend_mode = NI_EXTEND_NEAREST;
        }

        switch (extend_mode)
        {
        /* aaaaaaaa|abcd|dddddddd */
        case NI_EXTEND_NEAREST:
            src = first;
            dst = buffer;
            val = *src;
            while (size_before--)
            {
                *dst++ = val;
            }
            src = last - 1;
            dst = last;
            val = *src;
            while (size_after--)
            {
                *dst++ = val;
            }
            break;
        /* abcdabcd|abcd|abcdabcd */
        case NI_EXTEND_WRAP:
        case NI_EXTEND_GRID_WRAP:
            src = last - 1;
            dst = first - 1;
            while (size_before--)
            {
                *dst-- = *src--;
            }
            src = first;
            dst = last;
            while (size_after--)
            {
                *dst++ = *src++;
            }
            break;
        /* abcddcba|abcd|dcbaabcd */
        case NI_EXTEND_REFLECT:
            src = first;
            dst = first - 1;
            while (size_before && src < last)
            {
                *dst-- = *src++;
                --size_before;
            }
            src = last - 1;
            while (size_before--)
            {
                *dst-- = *src--;
            }
            src = last - 1;
            dst = last;
            while (size_after && src >= first)
            {
                *dst++ = *src--;
                --size_after;
            }
            src = first;
            while (size_after--)
            {
                *dst++ = *src++;
            }
            break;
        /* cbabcdcb|abcd|cbabcdcb */
        case NI_EXTEND_MIRROR:
            src = first + 1;
            dst = first - 1;
            while (size_before && src < last)
            {
                *dst-- = *src++;
                --size_before;
            }
            src = last - 2;
            while (size_before--)
            {
                *dst-- = *src--;
            }
            src = last - 2;
            dst = last;
            while (size_after && src >= first)
            {
                *dst++ = *src--;
                --size_after;
            }
            src = first + 1;
            while (size_after--)
            {
                *dst++ = *src++;
            }
            break;
        /* kkkkkkkk|abcd]kkkkkkkk */
        case NI_EXTEND_CONSTANT:
        case NI_EXTEND_GRID_CONSTANT:
            val = cval;
            dst = buffer;
            while (size_before--)
            {
                *dst++ = val;
            }
            dst = last;
            while (size_after--)
            {
                *dst++ = val;
            }
            break;
        default:
            break;
        }
    }

    Eigen::RowVectorXd output = Eigen::RowVectorXd::Zero(line_length);

    /* get lines: */
    double *iline = input_extend.data() + size1;
    double *oline = output.data();
    fw            = weights.data() + size1;
    /* the correlation calculation: */
    if (symmetric > 0)
    {
        for (int ll = 0; ll < line_length; ll++)
        {
            oline[ll] = iline[0] * fw[0];
            for (int jj = -size1; jj < 0; jj++)
                oline[ll] += (iline[jj] + iline[-jj]) * fw[jj];
            ++iline;
        }
    }
    else if (symmetric < 0)
    {
        for (int ll = 0; ll < line_length; ll++)
        {
            oline[ll] = iline[0] * fw[0];
            for (int jj = -size1; jj < 0; jj++)
                oline[ll] += (iline[jj] - iline[-jj]) * fw[jj];
            ++iline;
        }
    }
    else
    {
        for (int ll = 0; ll < line_length; ll++)
        {
            oline[ll] = iline[size2] * fw[size2];
            for (int jj = -size1; jj < size2; jj++)
                oline[ll] += iline[jj] * fw[jj];
            ++iline;
        }
    }

    return output;
}
/// @brief Calculate a 1-D convolution along the given axis.
/// The lines of the array along the given axis are convolved with the
/// given weights.
/// @param input ndarray
/// 1-D sequence of numbers.
/// @param weights ndarray
/// 1-D sequence of numbers.
/// @param mode
/// @param cval
/// @param origin
/// @return convolve1d : ndarray
/// Convolved array with same shape as input
Eigen::RowVectorXd convolve1d(Eigen::RowVectorXd input,
                              Eigen::RowVectorXd weights,
                              std::string mode = "reflect",
                              float cval       = 0.0,
                              int origin       = 0)
{
    // weights = weights[::-1]
    weights.reverseInPlace();
    origin = -origin;
    if (!(weights.size() & 1))
    {
        origin -= 1;
    }
    return correlate1d(input, weights, mode, cval, origin);
}
} // namespace ndimage

} // namespace scipy
